Split-merge cycle, fragmented collapse, and vortex disintegration in rotating 
Bose-Einstein condensates with attractive interactions 
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The dynamical instabilities and ensuing dynamics of singly- and doubly-quantized vortex states 
of Bose-Einstein condensates with attractive interactions are investigated using full 3D numerical 
simulations of the Gross-Pitaevskii equation. With increasing the strength of attractive interactions, 
a series of dynamical instabilities such as quadrupole, dipole, octupole, and monopole instabilities 
emerge. The most prominent instability depends on the strength of interactions, the geometry 
of the trapping potential, and deviations from the axisymmetry due to external perturbations. 
Singly-quantized vortices split into two clusters and subsequently undergo split-merge cycles in a 
pancake-shaped trap, whereas the split fragments immediately collapse in a spherical trap. Doubly- 
quantized vortices are always unstable to disintegration of the vortex core. If we suddenly change 
the strength of interaction to within a certain range, the vortex splits into three clusters, and one 
of the clusters collapses after a few split-merge cycles. The vortex split can be observed using a 
current experimental setup of the MIT group. 

PACS numbers: 03.75.Fi, 05.30.Jp, 32.80.Pj, 67.40. Vs 



INTRODUCTION 

Quantized vortices are among the hallmarks of super- 
fluidity, and have been observed in the gaseous Bose- 
Einstein condensates (BECs) of 87 Rb P,@ and 23 Na 0, 
where repulsive interatomic interactions play a crucial 
role in nucleation and stabilization of vortices. In the 
case of attractive interactions, however, creating vortices 
by stirring the condensate in a harmonic trap is diffi- 
cult because the vortex state is energetically unfavorable 
against decaying into the non-vortex ground state. A 
possible experimental scheme for creating a vortex state 
with attractive interactions is to prepare it in the repul- 
sive regime and then switch the sign of the interaction to 
attractive by using, e.g., the Feshbach resonance 
Another possible scheme to create vortices in an attrac- 
tive BEC is to imprint topological phases on the conden- 
sate using either controlled interconversion between two 
components Q, or spin rotation by adiabatic change 
of the magnetic field 0, • The aim of the present pa- 
per is to study the dynamical instabilities and ensuing 
dynamics of vortex states with attractive interactions in 
3D axisymmctric traps. 

There are two types of instabilities that cause a vortex 
state to decay: dynamical and thermodynamic instabil- 
ities. Dynamical instabilities, or modulational instabil- 
ities, set in when the Bogoliubov spectrum acquires an 
imaginary part which causes an exponential growth in 
the corresponding mode. Vortex nucleation |9l lld|. vor- 
tex exchange between two components ^1], splitting of 
a vortex with attractive interactions |l2l| , and disintegra- 
tion of a doubly-quantized vortex with repulsive interac- 
tions (l3lll4j are all attributed to dynamical instabilities. 
Thermodynamic instabilities are caused by interactions 



between the condensate and the thermal cloud accompa- 
nied by dissipation of the energy and angular momentum. 
Reductions in the number of vortices |2|, vortex bend- 
ing 0, 0] , and vortex-lattice formation are associ- 
ated with thermodynamic instabilities. The time scale 
of the latter instabilities at low temperatures is typically 
~ 1 s or longer In contrast, the time scale of 

dynamical instabilities can be much faster than that of 
thermodynamic instabilities since the former grows expo- 
nentially in time. We shall therefore focus on dynamical 
instabilities in this paper. 

In Ref. 0] , Pu et al. reported that singly- and doubly- 
quantized vortices in 2D exhibit a series of dynamical in- 
stabilities when the strength of interaction exceeds the 
corresponding critical values. We numerically solved the 
time-dependent Gross-Pitaevskii (GP) equation for a 2D 
system, and showed that the vortex splits into two clus- 
ters which revolve around the center of the trap, and un- 
dergo split-merge cycles or collapse [l^ . In the present 
paper, we examine such dynamics in full 3D systems 
and also study the dynamical instability of the doubly- 
quantized vortex, that was recently realized by the MIT 
group ||. Our primary findings are the following: 

(1) Our 3D analysis reveals that the dynamical insta- 
bilities depend crucially on the trap geometry. For a 
singly-quantized vortex, split-merge cycles can occur in 
a pancake-shaped trap, whereas the split fragments im- 
mediately collapse in a spherical trap. In a cigar-shaped 
trap, there is a parameter range in which the vortex can 
collapse into a single cluster due to the dipole instability. 

(2) The doubly-quantized vortex with attractive inter- 
actions is always unstable against disintegration of the 
vortex core. If the strength of interaction is suddenly 
changed to within a certain range, the vortex splits into 
three clusters which undergo split-merge cycles before 
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one or possibly more of them collapse. 

(3) Even when the lifetime of the condensate is very 
short, the phenomenon of vortex split can be observed 
by suddenly changing the interaction to a large attractive 
value. The split of a doubly-quantized vortex is shown to 
be observable using a current experimental setup of the 
MIT group. 

(4) Using a variational method, we find that both the 
split instability and the split-merge cycles can be under- 
stood by using only three basis functions. 

(5) In the immediate vicinity of the critical strength of 
interaction for the vortex split, quantum fluctuations like 
the ones found in a ID case jl9| become prominent. 

This paper is organized as follows. Section B briefly 
reviews the GP and Bogoliubov theories and discusses 
dynamical instabilities. Section D studies dynamical in- 
stabilities of a singly-quantized vortex state, and Sec. 
(hose of a doubly-quantized vortex state. Section B exam- 
ines a variational method to understand the dynamical 
instabilities and split-merge phenomenon. Section D shows 
that the split-merge phenomenon is accompanied by the 
emergence of a symmetry-broken low-lying state. Sec- 
tion B concludes this paper. 



FORMULATION OF THE PROBLEM 

We consider a dilute-gas BEC at zero temperature con- 
fined in an axisymmetric harmonic trap. The dynamics 
of the condensate are described by the GP equation 
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where m is the atomic mass, I4rap = mio\(x 2 + y 2 )/2 + 
mu> 2 /2 with u>± and lo z being the radial and axial trap 
frequencies, and a is the s-wave scattering length. Wc 
normalize the length, time, energy, and wave function 
by d = (h/mLU^) 1 / 2 , wj 1 , huj±, and (N/dl) 1 / 2 , respec- 
tively, where N is the number of atoms in the BEC, and 
the wave function is normalized as J dr\ip\ 2 = 1. The GP 
equation Q is then written as 

^^ = +\{r 2 + \ 2 z 2 H> + <?|VO, (2) 



where 



x 2 + y 2 , X — u z /uj_, and g = 4irNa/do. 



The Bogoliubov spectrum for a stationary state tpo 
of the GP equation is obtained from the Bogoliubov-de 



(K + 2g\iP \ 2 )v 
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where K = -V 2 /2 + (r 2 + \ 2 z 2 )/2 - ^ and ^ is the 
chemical potential. For the axisymmetric vortex states 
ipo ex e 1 ^, each Bogoliubov mode u n cx e l (i+ Sm )'t> couples 
only to v n oc e - l (i- Sm )4> ; where q and 5m are integers and 
<p = tan _1 y/a; is the azimuthal angle. In the following 
discussions, we shall refer to the modes with 8m = 1, 
2, and 3 as the dipole, quadrupole, and octupole modes, 
respectively. 

The stability of the system is characterized by the Bo- 
goliubov spectrum E n . When the system is in the ground 
state, or in a metastable state for which there is an en- 
ergy barrier that prevents the system from decaying into 
a lower energy state, all Bogoliubov eigenenergies are real 
and positive. An example of such a metastable state is 
a BEC with an attractive interaction below the critical 
strength of interaction for collapse [20I I2H Eif . When 
there are negative eigenvalues in the Bogoliubov spec- 
trum, the system can thermodynamically decay into a 
lower energy state by dissipating energy to the environ- 
ment. When complex eigenvalues emerge in the Bogoli- 
ubov spectrum, the corresponding modes grow exponen- 
tially in time and the system becomes dynamically unsta- 
ble. The time scale of such dynamical instability is given 
by the inverse of the imaginary part. The difference be- 
tween thermodynamic and dynamical instabilities is that 
the former requires energy dissipation, whereas the lat- 
ter can occur via an energy-conserving process. Thus, the 
dynamical instability can be a dominant decay process at 
sufficiently low temperature in a high-vacuum chamber 
in which dissipation is negligible over a short time. 

The time evolution of a system that deviates slightly 
from the stationary state xpo of the GP equation can be 
described in terms of the Bogoliubov modes u n and v n 
as [2 j 3j 



^ t +«e i ^ t ) 
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where £„'s (|e„| -€i 1) are small coefficients that describe 
the deviation. When the eigenvalue of a Bogoliubov 
mode becomes complex, the mode grows exponentially 
in time as 



J(q<t>-llt) 
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fa + ^ mEnt \e n \ [(«« + »„) cos(<5m0 - RcE n t + 6 n ) + {u n - v n ) sm{8m<P ~ ReE n t + 0„)] }, (5) 



where we leave only the Bogoliubov mode with the largest 
imaginary part ImE^ > 0, and define tpo = ipoe' lq ^, 



i(q-\-8m)4> 



v n e 
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Since both (u,v,E) and (u*,v*,E*) are solutions of 
Eq. ©, one of them satisfies lmE n > 0. Equation (J5J in- 
dicates that the <5m-fold density modulation grows expo- 
nentially in time while rotating at frequency HeE n /5m. 
Note, however, that the above discussion, which is based 
on Bogoliubov analysis, can be applied only to an initial 
stage of the growth of the unstable mode, and we must 
integrate GP equation (J2J to follow the time evolution 
over a longer time. 

We briefly describe the numerical procedure which 
will be used to solve Eqs. (0 and @. The stationary 
state ifjQ °f Eq. (0 is obtained by the Newton-Raphson 
method [13 or the imaginary time method [24| m the 
restricted functional subspace of ipo cx e lq ^ . The Bogoli- 
ubov spectra are then obtained by diagonalizing Eq. © 
in each subspace of the angular momentum 5m. Time 
evolution of the wave function is obtained by integrating 
Eq. (J2J using the Crank-Nicholson scheme [25j with a spa- 
tial discretization taken typically to be 256 x 256 x 256. 
In studying dynamical instabilities, we add a symmetry- 
breaking perturbation to the initial state such as 



tp(t = 0) oc (1 + 0.01 cosn</))^ , 



(6) 



where n is an integer. While dynamically unstable modes 
can grow from infinitesimal numerical noise over a suffi- 
ciently long time even when we start from an axisymmet- 
ric state, we add the perturbation JBJ in order to simulate 
experimental noise in a controllable manner so that the 
results are reproducible. This type of perturbation may 
be realized experimentally, e.g., using off-resonant laser 
beams propagating along the trap axis. The system is 
allowed to evolve in time until just before the collapse 
sets in, since numerical simulations of violent collapses 
and explosions jl^, in 3D space present a formidable 
computational challenge. 



DYNAMICAL INSTABILITIES IN A 
SINGLY-QUANTIZED VORTEX 

In this section, we study the dynamical instabilities 
in a singly-quantized vortex state ipo oc e*^. Figure Q 
shows imaginary parts of the Bogoliubov spectra as 
functions of the interaction parameter g for spherical 
(A = 1), pancake-shaped (A = 10), and cigar-shaped 
(A = 0.008 HI]) traps. As seen from Figs. □ (a) and 
(b), for the spherical and pancake-shaped traps, the 
quadrupole mode (8m = 2) first becomes unstable in the 
sense that the imaginary part ImE appears for values of 
\g\ smaller than those of other modes. A similar result is 
obtained in Ref. for a 2D system, which corresponds 
to the pancake-shaped trap with a large asymmetry pa- 
rameter A [29j |. For the cigar-shaped trap [Fig. ^ ( c )], 
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FIG. 1: Imaginary parts of the Bogoliubov spectra of the 
singly-quantized vortex state as functions of g for (a) A = 1 
(spherical trap), (b) A = 10 (pancake-shaped), and (c) A = 
0.008 (cigar-shaped). The dotted lines show critical values 
g"f below which the collapse of the vortex state is caused by 
the monopole instability. 



on the other hand, the dipole mode (5m = 1) first be- 
comes unstable, which occurs for A < 0.34 hj. This 
suggests that when the strength of interaction \g\ is adi- 
abatically increased, the dipole instability plays a domi- 
nant role of decay mechanisms for the cigar-shaped trap 
with A < 0.34; otherwise the quadrupole mode is domi- 
nant. 

Figurc|21shows the time evolution of integrated column 
densities and phase profiles for A = 10. The integrated 
column density is defined by 



Pi 



d £\ip\ 2 (£ = x,y,z), 



(7) 



and the phase of the integrated wave function J dzip is 
used to draw the phase profiles. The main panels in 
Fig. El show top views p z and the lower-left insets show 
side views p x . The interaction parameter is taken to 
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(d) t=18 |(e) t=25 1(f) t=14+l 




FIG. 2: (a)-(e) Snapshots of integrated densities p z (x,y) = 
J dz\ip\ 2 (main panels) and p x (y,z) — J dx\ip\ 2 (lower-left in- 
sets) for g — —6.5 and A = 10. The initial state is the sta- 
tionary singly-quantized vortex state for g = —6.5 disturbed 
by perturbation (JSJ with n — 2. The lower-right insets show 
gray-scale images of the phase of the integrated wave function 
<& z = J dzip, where the phase from (2n — l)n to (2n + l)ir is 
represented in gray scale from black to white, (f) The inte- 
grated density p z at t — 15, where the interaction is switched 
to g = 50 and the trap is turned off at t = 14. The sensitivity 
of the imaging in (f) is 12 times higher than that in (a)-(e). 
The sizes of the images are 6 x 6 in p z and the phase pro- 
files, and 6 x 2 in in (a)-(e), and 14 x 14 in (f) in units of 
(Ti/Wi) 172 - 



be g = —6.5, which exceeds the critical value for the 
quadrupole instability gn = —6.0, but is below the crit- 
ical value for the dipole instability g^ — —8.5. The Bo- 
goliubov eigenvalue of the quadrupole mode at g = —6.5 
is Eq ~ 1.5 + 0.40i. We add the perturbation JJjJ with 
n = 2 to the initial state to trigger the quadrupolar dy- 
namical instability. In Fig. [5] (b) , we see that a density 
fluctuation in the quadrupole mode starts to grow with 
the Lyapunov exponent IthEq ~ 0.4 and rotates at fre- 
quency KeEg/2 ~ 0.75 in accordance with Eq. (J5J. The 
vortex then splits into two clusters which revolve around 
the center of the trap [Fig. [5] (c)] , and subsequently these 
clusters unite to recover the original vortex [Fig. [21 (e)]. 
After this, the split-merge cycle repeats. Figure [21 (f) 
shows an expanded image of p z at t — 15, in which the 
interaction was switched from g = —6.5 to g = 50 and 
the trap was turned off at t = 14. The fringe pattern is 
due to interference between two matter waves emanating 
from the two clusters shown in Fig. [21 (c). The above 
results confirm that our findings in Ref. [12| , which were 
obtained in a 2D system, hold true in a pancake-shaped 
trap and are therefore realizable experimentally. 

When the vortex splits into n clusters, each of these 
subsequently collapses if the effective strength of inter- 
action per each cluster, i.e., \g\/n, exceeds the critical 
value Iffnonvortexl f° r the collapse of a condensate with- 
out vortices, while the split-merge cycle occurs if |<7|/n 




FIG. 3: Snapshots of integrated densities p z (x,y) (upper 
panels), p x (y,z) (lower panels), and phase profiles (insets) for 
g — —20 and A = 0.008. The initial state is the stationary 
singly-quantized vortex state for g — —20 disturbed by per- 
turbation (JSJ with n — 1. The sizes of the images are 6x6 
for p z and 6x4 for p x in units of (h/muj±) 1 '" 2 . 



is sufficiently smaller than |<7nonvortcxl- I n ^ ac ^ m the 
case of A = 10 in Fig. [21 \g\/2 — 3.25 is smaller than 
Iffnonvortexl = 3-73, and therefore the split clusters do not 
collapse, allowing the system to undergo split-merge cy- 
cles. However, for the spherical trap, ffnonvortcx = —7.22 
and (?q = —15.1, and hence the value |<?q|/2 = 7.55 at 
which the split occurs is always larger than IffnonvortexL 
indicating that split clusters always collapse. We per- 
formed numerical calculations for the spherical trap with 
\g\ slightly larger than \gg\, and confirmed that two clus- 
ters immediately collapse without undergoing a split- 
merge cycle (as shown in Figs. 4 (d)-(f) of Ref. |12j). 
The phenomenon, in which a vortex first splits into frag- 
ments and each fragment subsequently collapses, may be 
called "fragmented collapse" . 

Figure[31shows time evolution for the cigar-shaped trap 
with A = 0.008. The interaction parameter is taken to 
be g = —20, at which only the dipole mode is unsta- 
ble [see Fig. ^ (c)]. The initial state is the stationary 
state for g = — 20 plus the small perturbation 10 with 
n = 1. We see that the atoms tend towards one side 
of the ring due to the dipole instability [Fig. [21 (b)], af- 
ter which collapse occurs [Fig. (c)]. Since the dipole 
mode has angular momenta m = 1 ± 1 = 2 and 0, the 
phase defect that exists from the outset slightly devi- 
ates from the center associated with the growth of the 
to = component, and another phase defect enters due 
to the hi — 2 component as seen in the inset of Fig. [3 
(c). For A = 0.008, the critical value at which the dipole 
mode becomes unstable is g c ^ = —8.66, which exceeds 
fnonvortex = ~ 8.49, and therefore the dipole instability 
is always followed by collapse. We note that the dipole 
instability in Fig.^is distinct from the dissipative vortex 
spiral-out phenomenon that arises from the dipole mode 
with negative eigenvalue In the former, the 

quasiparticles that change angular momentum by ±1 are 
excited in pairs, and hence the angular momentum is 



conserved. The energy is also conserved, since the real 
parts of the excitation energies have the same magnitude 
with opposite signs. In the latter, the atoms are trans- 
ferred only into the non-vortex ground state decreasing 
the angular momentum and dissipating the energy. 



DYNAMICAL INSTABILITIES IN A 
DOUBLY-QUANTIZED VORTEX 

A doubly-quantized vortex state oc e 21 ^ of 23 Na BEC 
has recently been realized by the MIT group using a 
topological phase- imprinting method Q , where the inter- 
action is repulsive. The doubly-quantized vortex with re- 
pulsive interactions is predicted to be dynamically unsta- 
ble a gain st disintegration into two singly-quantized vor- 
tices |l3lll^l3^ |. In the case of attractive interactions, we 
will show that not only the vortex disintegration but also 
various new phenomena, such as a three-fold split-merge 
process, occur due to dynamical instabilities. 

The imaginary parts of the Bogoliubov eigenvalues for 
the doubly-quantized vortex states are shown in Fig. 0] 
We note that unlike the case of the singly-quantized vor- 
tex state, the quadrupole mode (dm — 2) is always dy- 
namically unstable for g < 0, indicating that no dynami- 
cally stable state exists for the case of attractive interac- 
tions. Therefore, an adiabatic decrease in g from positive 
to negative always gives rise to the quadrupole instability. 

Figure El shows time evolutions of the integrated den- 
sities p z (main panels) and p y (lower-left insets) and 
the phase profiles (lower-right insets) for the doubly- 
quantized vortex state in a spherical trap, where the 
initial state is a stationary state of the GP equation 
with g = — 20 plus the perturbation given in Eq. © 
with n — 2. For this strength of interaction, only the 
quadrupole mode is dynamically unstable [cf. Fig.|3](a)]. 
The vortex core disintegrates, as can be seen from the 
density and phase profiles. The angular momenta of the 
quadrupole excitations are m = 2 ± 2 = 4 and 0. Vortex 
disintegration is attributed to the growth of the m = 
component, and the four phase defects in Fig. [3(d) are 
due to the m — 4 component. The density distribution 
exhibits various changes such as high-density regions oc- 
curring in the opposite edge of the vortices [Fig. El (c)], 
then between the two vortices [Fig. 151(d)], and then the 
density distribution becomes similar to the split vortex 
as in the case of a singly-quantized vortex [Fig. E] (e)]. 
Finally, in Fig. [51 (f ) the central density exceeds critical 
density and collapses. 

Another important distinction of Fig. 0] from Fig. 
is that there is a range of g over which the imaginary 
part of the eigenvalue of the octupole mode (dm = 3) be- 
comes largest. This implies that if g is suddenly changed 
into that regime, the octupolar instability dominates un- 
til the quadrupole mode grows. Figure El shows a typical 
example of such situations with A = 10, where the initial 
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FIG. 4: Imaginary parts of the Bogoliubov spectra of the 
doubly-quantized vortex state as functions of g for (a) A = 1, 
(b) A = 10, and (c) A = 0.008. The dotted lines show the 
critical value g^j below which the collapse of the vortex state 
is caused by the monopole instability. 



state is the noninteracting stationary state with vorticity 
of two subject to perturbation JSJ with n = 3, and the 
interaction parameter is changed from g = to g = — 10 
at t = 0. As the octupole mode grows, the vortex splits 
into three clusters which revolve around the center of the 
trap [Fig. El (c)]. Interestingly, the three clusters merge 
to recover the original vortex [Fig. (d)] as in the two- 
cluster case shown in Fig.|2 This is because the effective 
strength of interaction per each cluster |<?|/3 = 3.33 is 
smaller than the critical value for collapse of the non- 
vortex state Iffnonvortexl = 3.73. Subsequently, the vor- 
tex once again splits into three clusters, and the balance 
in atomic numbers between these is then broken by the 
quadrupole instability. It follows that one cluster [de- 
noted by A in Fig. (f)] grows, then another cluster [B 
in Fig. El (g)] grows, and finally the remaining cluster [C 
in Fig. El (h)] grows and eventually collapses. This phe- 
nomenon is similar to the seesaw-like oscillation between 
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FIG. 5: Snapshots of integrated densities p z (x, y) (main pan- 
els), p y (x, z) (lower-left insets), and phase profiles (lower-right 
insets) for A = 1. The initial state is the stationary doubly- 
quantized vortex state for g — disturbed by perturbation 
© with n — 2, and the interaction is switched to g = —20 at 
t = 0. The sizes of the images are 6x6 for p z and the phase 
profiles, and 6 x 3 for p y in units of (Ti/mLo^) 1 ^ 2 . 



two clusters due to the dipole instability |l2( . At the ini- 
tial stage of disintegration in which the octupole mode 
plays a dominant role [Figs. (a) -(e)], the split in the 
phase defects is not very prominent because the angular 
momenta of the Sm = ±3 modes are m = 2 ± 3 = 5 and 
— 1. When the quadrupole mode becomes significant, the 
outward shifts of the phase defects manifest themselves 
[Figs. (g)-(h)] due to the m — component of the 
quadrupole mode. FigureEl(i) shows the integrated den- 
sity at t — 13, following the interaction being switched 
from g = — 10 to g = 50 and the trap being turned off 
at t = 12. We see the three- fold fringe pattern resulting 
from interference between the matter waves that emanate 
from the three clusters shown in Fig.El(c). 

In an elongated cigar-shaped trap, dynamics along the 
trap axis also plays an important role. Figure El shows 
the integrated density and phase profiles for A = 0.008. 
The initial state is the noninteracting stationary state 
with vorticity of two plus perturbation JSJ) with n = 2, 
and the interaction is switched from g = to g = —38 
at t = 0. The atoms first gather along the symmetry 
axis (i.e., the z axis) towards z = [Fig. Etb)], then 
the vortex begins to disintegrate since the density is high 
enough. Subsequently the atoms rebound along the z di- 
rection [Fig. Etc)]. We see that the disintegrated vortex 
lines in Fig. El (c) are straight, unlike the repulsive case 
where twisting of vortex lines occurs 0|. In Fig. El (d), 
the condensate splits into two high-density regions con- 
taining the disintegrated vortices as a consequence of the 
rebound. 

In an experiment using 23 Na, the lifetime of the con- 
densate is as short as 1 or 2 ms near the Feshbach res- 
onance [iiij], which indicates that the lifetime is t < 3 
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FIG. 6: (a)-(h) Snapshots of integrated densities p z {x,y) 
(main panels), p x {y,z) (lower- left insets), and phase profiles 
(lower-right insets) for the pancake-shaped trap with A = 10. 
The initial state is the stationary doubly-quantized vortex 
state for g = disturbed by perturbation J§J with n = 3, 
and the interaction is switched to g — — 10 at t = 0. (i) 
The integrated density p z at t — 13, where the interaction 
is switched to g — 50 and the trap is turned off at t = 12. 
The sensitivity of the imaging is 12 times higher than that in 
(a)-(h). The sizes of the images are 6 x 6 in p z and the phase 
profiles, and 6 x 2 in p x in (a)-(h), and 14 x 14 in (i) in units 
of (h/m,Lo ± ) 1/2 . 



(in our dimensionless units of time) for the optical dipole 
trap used in the MIT grou p (1 1 = 27r x 250 Hz, lu z — 
2ir x 2 Hz, and A = 0.008) [23]. For dynamical instabil- 
ities to occur within this short lifetime, the strength of 
attractive interaction \g\ must be much larger than those 
in the above cases. Figure |S] shows the density and phase 
profiles just before the collapse, where the initial state 
is the noninteracting stationary state with vorticity of 
two plus perturbation |JBJ with n = 1, 2, 3, and 4, and 
the interaction is switched from g = to g = —500 at 
t = 0. Wc find that the collapsing dynamics depend on 
the initial perturbation n, since the system has several 
multipole instabilities for g = —500. We note that the 
vortex split occurs only around z — 0, where the density 
is high enough. Thus, the vortex-split phenomena can be 
observed for large values of \g\, even when the lifetime of 
the condensate is very short. 
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FIG. 7: Snapshots of p z (x,y) (upper panels), phase profiles 
(insets), p x (y,z) (lower-left panels), and p y (x,z) (lower-right 
panels) for the cigar-shaped trap with A = 0.008. The initial 
state is the stationary doubly-quantized vortex state for g = 
disturbed by perturbation © with n — 2, and the interaction 
is switched to g — —38 at t — 0. The sizes of the images are 
6 x 6 for p z and the phase profiles, and 6 x 16 for p x and p y 
in units of (h/mu;±) 1 ' 12 . 



VARIATIONAL ANALYSIS 

To understand the split-merge phenomena shown in 
the preceding sections, we restrict ourselves to 2D space 
in this and next sections. 

We employ a Gaussian variational wave function for 
the condensate as 



^q 



M 



r+iq 



+i 



(8) 



where q is the vorticity and d is a variational parameter 
characterizing the size of the condensate. Substituting 
Eq. JSJ) into the GP energy functional 



EM 



dr 



and minimizing it with respect to d gives 

2%l+l|g|!2(l + | g |) 7r 



1 



(9) 



(10) 



For g < gfj = -2 2 l'l +1 |g|! 2 (l + |g|)7r/(2|g|)!, the energy 
functional Q has no extremum, and the system collapses 
due to the monopole instability. The critical values g^ = 
— 8n for q = 1 and = — 16n for q = 2 agree reasonably 
well with the numerically obtained critical values —24 
and —45, respectively. 

To study the time evolution from Eq. JSJ, we use a 
variational wave function given by 



(n; 



where 



1 and ifj m is given by Eq. © with 



q replaced by m. In Eq. I|ll[) we use the same width 



(a) (b) (c) (d) 
n=l n=2 n=3 n=4 
t=2.7 t=2.2 i=1.9 i=22 





FIG. 8: Snapshots of p z (x,y), phase profiles, and p y (x,z) 
(from top to bottom) for the cigar-shaped trap with A = 
0.008. The initial state is the stationary doubly-quantized 
vortex state for g — disturbed by perturbation © with (a) 
n = 1, (b) n = 2, (c) n = 3, and (d) n — 4. The interaction 
is switched from g — to g — —500 at t = 0. The sizes of the 
images are 6x6 for p z and the phase profiles, and 6 x 20 for 
p y in units of {Ti/muj± 



U/2 



d as Eq. l(TU|) for all -0m and neglect the radial degrees 
of freedom to simplify the following calculation and to 
obtain the results analytically. 

Substituting Eq. I|llfl into the action 



d 

dtE[ip] - i I drdtip*—ip 



(12) 



yields 



S = 



dt —i > c m c 



Er<m\, 
^3, 



7712 * * 

m 4 c m 1 C m 2 Cm 3 C ™4 



(13) 



77117712 
77I377I4 



7711,7712 — 
7713,7714 



where e m = (1 + |m|)(d 2 + l/d 2 )/2, and G 
/ dxi\)* mi ^* m ^l) m ^ mi , which vanishes for mi+m 2 -m 3 - 
m 4 7^ 0. From dSjdc* m = 0, the variational equation of 
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motion for c m becomes 



c r J- n \ rim,m-i * 



(14) 



mim2rri3 



We assume |c g | ~ 1 and \c m ^ q \ <C 1, and neglect terms of 
the order of 0(\c m -^ q \ 2 ). Equation (fH|) for m = q reduces 
to 



(ex + 



(%|)!fl 



2 2 l9l+ 1 | 9 |! 2 7rd 2 ' 



(15) 



and hence we obtain 



defining c„ 
tions 



a e~^ lt . Using Eq. and 
= e -t,i *Cm, we obtain a closed set of equa- 



/c ,, 



-29 -m 



(e ro - m + 2 5 G™;«) c m 

J-nf7 m ' 2<i, ~ m r* 

(s2q-m — M + ^gG^l-m l) <S 



(16a) 



2q—m 



y q.q 4ii 



(16b) 



from which the eigenfrequencies can be found. The con- 
dition under which the eigenfrequencies are real is given 
by 



A, 



q,2q—m 



e m + e 2q - m -2» + 2g (g^ + G^™ 
? 2 (G™ g 2? -" 1 ) 2 > 0. 



(17) 



First we consider the singly-quantized vortex (q = 1) 
Stability of the quadrupole mode is determined by 



Di, 



5g 
8tt 



1 



1 + JL - 



1 



9_ 

2tt 



32tt ; 



(18) 



which becomes negative for g < —9.2. This value is in 
reasonable agreement with the numerically obtained crit- 
ical value for the quadrupole mode in 2D, Qq = —7.8. 
However, the discriminant D\$ — g 2 / (64ir 2 d 4 ) is al- 
ways positive, which contradicts the numerical result that 
dipole instability does occur for g < gg = —11.5. Also, 
for m = —2 or m = 4, the numerical calculation shows 
that dynamical instability occurs for g < —16.9, while 
D\ t & is always positive for g > —Sir [g must be larger than 
— 8n as seen from Eq. (fTTTf) ] . In the case of the doubly- 
quantized vortex (q = 2), D 2 , 4 = -475 2 /(4096tt 2 g? 4 ) is 
always negative, which is consistent with the fact that 
the quadrupole mode is always dynamically unstable for 
g < 0. However, Z3 2 ,4 fails to reproduce the behavior for 
g > 0, where stable and unstable regions appear alterna- 
tively 0,0]. The discriminant D 2i5 becomes negative 
for g < —15.8, which is in reasonable agreement with 
the numerical result that the octupole instability sets in 
for g < -11.7. However, D 2 ,3 = g 2 / '(256tt 2 d 2 ) > con- 
tradicts the numerical result that dipole instability does 
occur for g < g c £ = —19.8. 




FIG. 9: Time evolutions of Ci (solid curve), C3 (dashed curve), 
and c_i (dotted curve) in the three- mode approximation of 
Eq. H14H for g — —12. The initial condition is C3 = c_i = 0.01 
and ci = (1 — 2 • O.Ol 2 ) 1 ^ 2 . The insets show the density and 
phase profiles at t = and t = 7.7, where the sizes of the 
images are 6 x 6 in units of (h/muj) 1 ^ 2 . 



Thus, the Gaussian ansatz Ijllfl is viable for some cases, 
whereas more sophisticated trial functions are needed to 
completely describe the dynamical instabilities. The vari- 
ational results can be improved by using the different 
variational width d for each ip m in Eq. IjllJI , which may be 
determined by, e.g., the method used in Ref. j^. How- 
ever, near and above the critical point, this variational 
method fails, since the norm of the Bogoliubov function 
vanishes when a complex eigenvalue emerges. 

Equation l|18fl is expanded near the critical point as 



Di, 3 = 



3V3 



(4%/3-3V2)7r 



(.9 ~ .9") + O {(g - <? CI ') 2 ) , (19) 



where g cr = 16(\/6 — 3)tt/3. Just above the critical point, 
therefore, the imaginary part of the complex eigenvalue 
is proportional to (-L>i, 3 ) 1/2 oc {g CI - g) 1 / 2 [l^. All the 
imaginary parts shown in Figs. ^ and 01 (except 8m = 2 
in Fig.^J have this g dependence near the critical points 
[so do the 6m = 1 curves in Figs. ^ ( c ) and 01 (c) in the 
immediate vicinity of the critical points]. The 6m = 2 
lines for q = 2 in Fig.^Jare proportional to \g\ near g = 0, 
which agrees with (— -D^) 1 / 2 oc \g\. 

The split-merge cycles shown in Figs. |5| and are re- 
currence phenomena arising from nonlinearity. To under- 
stand the behavior in Fig. El qualitatively, we reduce the 
number of modes in Eq. (|14(l to three, i.e., m = 1, 3, and 
— 1. Equation (|14|) then reduces to simultaneous nonlin- 
ear differential equations with three variables, which can 
be solved numerically. Figure shows the time evolution 
of |ci| 2 , | C3 1 2 , and |c_i| 2 for g = —12, where the initial 
condition is c 3 = c_i = 0.01 and c x = (1 - 2 • 0.01 2 ) 1/2 . 
The vortex can be seen to split at t ~ 7.7, but subse- 
quently recovers the original vortex during 12 < t < 18, 
and then splits again at t ~ 23.2. Thus, the split-merge 
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phenomenon can be qualitatively reproduced by a simple 
three-mode approximation. 



LOW-LYING SOLITON STATE AND 
GOLDSTONE MODE 



The split-merge dynamics shown in Fig. [21 implies the 
existence of a low-lying state in which two clusters revolve 
around each other without changing their shapes. To find 
such a low-lying stationary state, we numerically seek 
a local minimum of the GP energy functional with the 
angular momentum held fixed at that of the initial vortex 
state (L z ) = 1, where L z = —id/d<fi. The GP energy 
functional with Lagrange multipliers is given by 



(a) 2 



F[»p] = E[i/j] + / dr 



(20) 



where E[ip] is given in Eq. 10, and the chemical potential 
/i and the angular frequency fl of rotation of the system 
are introduced as the Lagrange multipliers in order to fix 
the normalization and the angular momentum. 

Using the Newton-Raphson method starting from 
appropriate initial states, we obtain the stationary states 
in the frame of reference rotating at frequency Q. Fig- 
ure 1101 shows various properties of the axisymmetric 
vortex state and the low-lying state. When \g\ is in- 
creased, the stationary state bifurcates at the critical 
point <7q = —7.79 of the quadrupole instability, at which 
the low-lying state emerges as represented by the lower 
branches of E and fi in Fig. EH (a). Density profiles of 
the low-lying states for g = — 8 and g = —10 are shown 
as gray-scale images in Fig. EH ( a )- It is clear from them 
that the split-merge cycle is due to oscillations between 
the axisymmetric vortex state (upper branch) and the 
low- lying soliton state (lower branch) . It is interesting to 
note that the energies of the axisymmetric and low-lying 
states are smoothly connected at the critical point, while 
the chemical potential has a kink, which implies that the 
nature of the phase transition between the vortex state 
and the soliton state is of second order. 

The excitation energies are plotted in Fig. ^3 (b), 
where E$ m denotes the energy of the lowest excitation 
with angular momentum 1 + 5m above the axisymmetric 
vortex state. As g approaches Qq = —7.79 in the axisym- 
metric regime, the excitation energy Es m =2 + Es m =-2 
vanishes, despite the fact that both Es m =2 and Es m =2 
are nonzero at <?q. This implies that quasiparticles with 
dm = ±2 can be excited in pairs without changing en- 
ergy and angular momentum. Hence, the two modes 
dm = ±2 may be strongly coupled via the pair exci- 
tations. If we follow the axisymmetric branch [the up- 
per branch in Fig. ( a )] , the eigenenergies of the two 
modes become complex. If we follow the low-lying soli- 
ton branch, one of the two modes becomes the Goldstone 



a. 




(b) 2 
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FIG. 10: (a) The energy E (solid curves) and chemical poten- 
tial /i (dashed curves) of the axisymmetric singly-quantized 
vortex state (upper branch) and low-lying soliton state (lower 
branch for g < —7.79). Gray-scale images show the density 
profile of the axisymmetric vortex state for g — —7.7 and 
those of the low-lying states for g = —8 and g — —10. (b) 
The Bogoliubov excitation energies of the quadrupole modes 
above the axisymmetric vortex state (g > —7.79) and those 
of the corresponding modes above the low-lying soliton state 
(g < —7.79). (c) The number of virtually excited atoms in the 
Bogoliubov modes. The inset is the logarithmic plot of the 
excitation energies and the number of excited atoms (both of 
them refer to the common scale) as functions of \g — Qq\. The 
vertical dotted lines in (a)-(c) show gg = —7.79. 



mode associated with the axisymmetry breaking, and the 
energy of the other mode E\ increases with decreasing g. 

Figure 1101 (c) shows the number of virtually excited 
atoms J dr\v\ 2 in the quadrupole modes Ns m =2+Ns m =-2 
for g > <7q and that in the corresponding mode Aq for 
g < g C Q. At g*Q the number of excited atoms diverges, 
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indicating that the Bogoliubov approximation breaks 
down. The inset of Fig. ^| (c) shows the logarithmic 
plot of the excitation energies and the numbers of excited 
atoms. Close to the critical point \g — <7g| <C 1, the exci- 
tation energies are found to be proportional to \g — ffgl 1 ^ 2 
and the numbers of excited atoms to \g — .9q| -1 ^ 2 . 

We note that these behaviors near the critical point 
(the kink in /i, the softening of the excitation mode, the 
divergence in the number of excited atoms, and their 
power laws) are very similar to those occurring in the 
transition between the uniform-density and soliton states 
in a ID ring [lflj . Thus, in the immediate vicinity of 
the critical point at which the vortex split occurs, the 
mean-field and Bogoliubov approximations break down 
due to large quantum fluctuations, and we must em- 
ploy, e.g., the exact diagonalization method 0] to study 
this regime. However, the mean-field theory used in the 
present paper is still valid except for this narrow critical 
region. 



CONCLUSIONS 

We studied dynamical instabilities in singly- and 
doubly-quantized vortex states with attractive interac- 
tions, and found a rich variety of dynamics triggered by 
the dynamical instabilities. The dynamics depend on the 
strength of interaction g, the trap geometry A, the initial 
perturbations, and the manner (adiabatic or sudden) in 
which g changes. 

A singly-quantized vortex in a trap with A > 0.34 first 
exhibits the quadrupole instability as g is adiabatically 
decreased. The vortex then splits into two clusters re- 
volving around the center of the trap. In a spherical trap 
the split fragments immediately collapse for any value 
of g that gives rise to the split instability, whereas in a 
pancake-shaped trap with A = 10 they unite to recover 
the original vortex and then the split-merge cycles repeat 
for \g\ slightly exceeding \gq\- In a cigar-shaped trap 
with A < 0.34, the dipole instability causes the atoms to 
concentrate in one cluster and then the system collapses. 

The doubly-quantized vortex state is always dynami- 
cally unstable against disintegration of the vortex core 
for g < 0, and therefore the disintegration is unavoidable 
if we adiabatically decrease g. If g is suddenly changed to 
some negative value with an appropriate initial perturba- 
tion to the condensate, the octupole instability can domi- 
nate the quadrupole instability, and the vortex splits into 
three clusters. In a pancake-shaped trap, these undergo 
split-merge cycles as in the case of a singly-quantized 
vortex, but eventually one of them collapses due to the 
quadrupole instability that always exists. We have also 
demonstrated that the split phenomena can be observed 
using the current experimental setup of the MIT group. 

To understand the dynamical instability and the split- 



merge cycles, we performed a variational analysis. The 
emergence of a complex eigenvalue in some modes was 
reproduced by the Gaussian trial functions. We were 
able to reproduce the split-merge cycles by assuming only 
three relevant modes with angular momenta 1 and 1 ± 2, 
suggesting that only a few modes are associated with the 
split-merge cycles. 

In the immediate vicinity of the critical strength of in- 
teraction of the dynamical instability, the transition from 
the axisymmetric vortex state to the split state is accom- 
panied by divergence in the number of virtually excited 
atoms and by softening of the excitation mode according 
to the Bogoliubov theory, which is very similar to the ID 
case discussed in Ref. |l9j . The dynamical instabilities in 
vortex states of attractive BECs therefore involve a fun- 
damental generic issue of the decay of a many-particle 
quantum system. 

Vortex-split phenomena are attributed to the interplay 
between the attractive interaction and the topological 
phase defect, that is, the condensate is bound to break 
axisymmetry, since they cannot gather at the phase de- 
fect. The emergence of the dynamical instabilities and 
ensuing dynamics may be regarded as pattern formation 
due to symmetry breaking. More complicated and in- 
teresting pattern might be formed in the vortex lattice, 
which will be published elsewhere. 

Dynamics after the collapse were not studied because 
3D simulation of the collapse and explosion phenomena 
exceeds our present computational ability. The burst 
atoms ejected from each cluster interfere with each other 
and exhibit some anisotropic pattern if the burst atomic 
cloud is coherent. Therefore experiments on vortex split 
and fragmented collapse will provide a test of whether or 
not the burst atoms are coherent. 

We thank J. K. Chin for discussions. This work was 
supported by the Special Coordination Funds for Pro- 
moting Science and Technology and a Grant-in-Aid for 
Scientific Research (Grant No. 15340129) by the Min- 
istry of Education, Science, Sports, and Culture of Japan, 
and by the Yamada Science Foundation. 
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